#Figure 3: Ratio %VA-Backed Loans Granted to White Applicants to %White Veterans in 25 Metropolitan Areas (1950)
library(tidyverse)
library(readxl)
library(here)
library(sf)  
library(tigris)
options(tigris_use_cache = TRUE) 
library(haven)
library(ggplot2)
library(usmap)
library(ggspatial)
library(ggmap)
library(maps)
library(mapdata)
library(dplyr)
library(readxl) 
library("stringr")
library(foreign)


setwd ("YourFilePath")
getwd()

# read in CBSA's(core-based statistical areas ) from census
usmetro<-core_based_statistical_areas( year = 2018) 
list(usmetro$GEOID, usmetro$NAME) 
usmetro$state <- str_sub(usmetro$NAME,-2,-1) 

# restrict the cbsas to contiguous states 
contig_usmetro <- subset(usmetro, state !="AK" & state != "HI" & state !="PR") 
class(contig_usmetro) #"sp"

#Read in us state outlines
usstate<- states()
list(usstate$NAME) #
usstate<- subset(usstate, NAME != "Commonwealth of the Northern Mariana Islands" & NAME != "United States Virgin Islands"  & NAME != "American Samoa" & NAME !=  "Guam" & NAME != "Puerto Rico"  & NAME != "Hawaii" & NAME != "Alaska")

#change us state outlines into sf class for ggplot 
usstate <- usstate %>%  
  st_as_sf(coords = c("INTPTLON", "INTPTLAT")) 
class(usstate) 
st_geometry(usstate) 
class(usstate)

list(usstate$geometry)

# read in the va loans data
valoans <- read_excel("map_HLG.xls") 
valoans <- select(valoans, ratio_loanpop, ratio_whitevavetii, NAME) 
class(valoans) #tbl

#merge loans with all the metro area in us (contig_usmetro)
valoans_cbsas <- merge(contig_usmetro, valoans, by = "NAME") 
summary(valoans_cbsas$ratio_whitevavetii) 
names(valoans_cbsas) 
class(valoans_cbsas)

# Transform valoans_cbsas to sf
valoans_cbsas <- valoans_cbsas %>%  
  st_as_sf(coords = c("INTPTLON", "INTPTLAT"))
class(valoans_cbsas) 
st_geometry(valoans_cbsas) 

# Figure 3: Ratio %VA-Backed Loans Granted to White Applicants to %White Veterans in 25 Metropolitan Areas (1950)
valoans_cbsas$ratio_vetbins <- cut(valoans_cbsas$ratio_whitevavetii, c(0.95, 1.00, 1.05, 1.10, 1.15, 1.20, 1.25))
select(valoans_cbsas, ratio_vetbins, ratio_whitevavetii)
summary(valoans_cbsas$ratio_vetbins) 

ggplot() +
  geom_sf(data = usstate, size=0.5, color = "black", fill = "transparent") +
  geom_sf(data = contig_usmetro, size = 0.2, color = "black", fill = "transparent") +
  geom_sf(data = valoans_cbsas, size=0.15, color = "black", aes(fill = as.factor(ratio_vetbins))) + 
  scale_fill_brewer(palette="YlGnBu", na.value="white", 
                    labels = c("0.95 - 0.99", "1.00 - 1.05", "1.06 - 1.10", "1.11 - 1.15", "1.16 - 1.20", "1.21 - 1.25")) +  
  theme(
    panel.border = element_blank(),     
    panel.grid.major = element_blank(),    
    panel.grid.minor = element_blank(),   
    panel.background = element_blank(), 
  ) 

ggsave("fig3_map.png") 




